clear

cd ../china/data
load hs10_indlist.raw;
indN=length(hs10_indlist);

Mean_ProbChange_priceshock_all=zeros(1,1)
Mean_ProbChange_priceshock_all_D2=zeros(1,1)
Mean_ProbChange_qualshock_all=zeros(1,1)
Mean_ProbChange_qualshock_all_D2=zeros(1,1)

sample_select=zeros(51,1);
for y=1:50
if (y>1 & y<5)
sample_select(y,1)=1;
else
if (y>5 & y<43)
sample_select(y,1)=1;
end
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Start by calculating Average Marginal Effects %%%%
%% d P(stay) / d p  and d P(stay / d lambda %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for w=1:50

if sample_select(w,1)==1
%if (w>1 & w<43)

folder=sprintf('../china/estimation/hs10/regular/raw_%d',hs10_indlist(w,1))


cd (folder)


load imp.raw;
load data2_q_regular.mat
load beta_regular.mat;
load EV.mat

else

folder=sprintf('../china/estimation/hs10/biggest/raw_%d',hs10_indlist(w,1))


cd (folder)


load imp.raw;
load data2_q_biggest.mat
load beta_biggest.mat;
load EV.mat

end

load imp.raw;
load price_minyear.raw;
load cities.raw

price_sd=sqrt(var(price_minyear));
lambda_sd=sqrt(var(cities(:,2)));

nExp = X;
nImp = M;
S=N*nExp*nExp;
delta=0.975;

p=(1:N)';
p=repmat(p,X*X,1);

x=1:X;
x=repmat(x,N,1);
x=x(:);
x=repmat(x,X,1);

xm1=1:X;
xm1=repmat(xm1,N*X,1);
xm1=xm1(:);

[beta_p_sol beta_x_sol beta_c_sol xi_tilde_sol];


U=beta_p_sol*ep + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*lambda;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_0 = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;


% Everyone stays with their partner
impprice=imp(:,1);
impxm1=imp(:,2);
impxm0=imp(:,3);
impxm0=impxm1;

statenum=[1:S]';

stay_state=zeros(M,1);

for h=1:M
    for i=1:N
        if impprice(h,1)==i
            for j=1:X
                if impxm0(h,1)==j
                    for k=1:X
                        if impxm1(h,1)==k
                            stay_state(h,1)=statenum((k-1)*X*N+(j-1)*N+i,1);
                        end
                    end
                end
            end
        end
    end
end

P_0_stay = P_0 ( stay_state);

% Shock price at original supplier

new_ep = ep;
for h=1:S
  if xm1(h,1) == x(h,1)
      new_ep(h,1) = ep(h,1) + price_sd;
  end
end


U=beta_p_sol*(new_ep) + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*lambda;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_priceshock = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;

P_priceshock_stay = P_priceshock ( stay_state);


Mean_ProbChange_priceshock= mean((P_priceshock_stay-P_0_stay))

[P_0_stay P_priceshock_stay];

Mean_ProbChange_priceshock_all(w,1)=Mean_ProbChange_priceshock;

% Shock quality at original supplier

new_lambda = lambda;
for h=1:S
  if xm1(h,1) == x(h,1)
      new_lambda(h,1) = lambda(h,1) + lambda_sd;
  end
end


U=beta_p_sol*ep + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*(new_lambda);

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_qualshock = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;

P_qualshock_stay = P_qualshock ( stay_state);


Mean_ProbChange_qualshock= mean((P_qualshock_stay-P_0_stay))

[P_0_stay P_qualshock_stay];

Mean_ProbChange_qualshock_all(w,1)=Mean_ProbChange_qualshock;

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate SEs for Average Marginal Effects using Delta Method %%%%
%% For this we need d^2 P(stay) / d beta_p d p 
%% and d^2 P(stay) / d xi d lambda
%% Then take sqrt of those things times var(beta) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd ../china/estimation/standard_errors
load standard_errors

for w=1:50

if sample_select(w,1)==1
%if (w>1 & w<43)

folder=sprintf('../china/estimation/hs10/regular/raw_%d',hs10_indlist(w,1))


cd (folder)


load imp.raw;
load data2_q_regular.mat
load beta_newq_regular.mat;
load EV_newq.mat

else

folder=sprintf('../china/estimation/hs10/biggest/raw_%d',hs10_indlist(w,1))


cd (folder)


load imp.raw;
load data2_q_biggest.mat
load beta_newq_biggest.mat;
load EV_newq.mat

end

load imp.raw;
load price_minyear.raw;
load cities.raw

price_sd=sqrt(var(price_minyear));
lambda_sd=sqrt(var(cities(:,2)));

nExp = X;
nImp = M;
S=N*nExp*nExp;
delta=0.975;

p=(1:N)';
p=repmat(p,X*X,1);

x=1:X;
x=repmat(x,N,1);
x=x(:);
x=repmat(x,X,1);

xm1=1:X;
xm1=repmat(xm1,N*X,1);
xm1=xm1(:);

[beta_p_sol beta_x_sol beta_c_sol xi_tilde_sol];

%%%%% Re-do with beta_p shocked too

new_beta_p = beta_p_sol + standard_errors(w,1)


U= new_beta_p*ep + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*lambda;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_0 = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;


% Everyone stays with their partner
impprice=imp(:,1);
impxm1=imp(:,2);
impxm0=imp(:,3);
impxm0=impxm1;

statenum=[1:S]';

stay_state=zeros(M,1);

for h=1:M
    for i=1:N
        if impprice(h,1)==i
            for j=1:X
                if impxm0(h,1)==j
                    for k=1:X
                        if impxm1(h,1)==k
                            stay_state(h,1)=statenum((k-1)*X*N+(j-1)*N+i,1);
                        end
                    end
                end
            end
        end
    end
end

P_0_stay = P_0 ( stay_state);

% Shock price at original supplier

new_ep = ep;
for h=1:S
  if xm1(h,1) == x(h,1)
      new_ep(h,1) = ep(h,1) + price_sd;
  end
end


U=new_beta_p*(new_ep) + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*lambda;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_priceshock = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;

P_priceshock_stay = P_priceshock ( stay_state);


Mean_ProbChange_priceshock_D2= mean((P_priceshock_stay-P_0_stay))

[P_0_stay P_priceshock_stay];

Mean_ProbChange_priceshock_all_D2(w,1)=Mean_ProbChange_priceshock_D2;

%%%%%%%%% Re-do with xi shocked too.

new_xi = xi_tilde_sol+ standard_errors(w,4);

U=(beta_p_sol)*ep + beta_x_sol*switched + beta_c_sol*city_switched + (new_xi)*lambda;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_0 = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;


% Everyone stays with their partner
impprice=imp(:,1);
impxm1=imp(:,2);
impxm0=imp(:,3);
impxm0=impxm1;

statenum=[1:S]';

stay_state=zeros(M,1);

for h=1:M
    for i=1:N
        if impprice(h,1)==i
            for j=1:X
                if impxm0(h,1)==j
                    for k=1:X
                        if impxm1(h,1)==k
                            stay_state(h,1)=statenum((k-1)*X*N+(j-1)*N+i,1);
                        end
                    end
                end
            end
        end
    end
end

P_0_stay = P_0 ( stay_state);



% Shock quality at original supplier

new_lambda = lambda;
for h=1:S
  if xm1(h,1) == x(h,1)
      new_lambda(h,1) = lambda(h,1) + lambda_sd;
  end
end


U=beta_p_sol*ep + beta_x_sol*switched + beta_c_sol*city_switched + (new_xi)*(new_lambda);

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



P_qualshock = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;

P_qualshock_stay = P_qualshock ( stay_state);


Mean_ProbChange_qualshock_D2= mean((P_qualshock_stay-P_0_stay))

[P_0_stay P_qualshock_stay];

Mean_ProbChange_qualshock_all_D2(w,1)=Mean_ProbChange_qualshock_D2;


end

Av_Marg_Effect_price=zeros(1,1)
Av_Marg_Effect_price_SE=zeros(1,1)
Av_Marg_Effect_qual =zeros(1,1)
Av_Marg_Effect_qual_SE =zeros(1,1)

for w=1:51
Av_Marg_Effect_price(w,1)    = Mean_ProbChange_priceshock_all(w,1)
Av_Marg_Effect_price_SE(w,1) = sqrt( Mean_ProbChange_priceshock_all_D2(w,1) * standard_errors(w,1)^2 * Mean_ProbChange_priceshock_all_D2(w,1) )
Av_Marg_Effect_qual(w,1)     = Mean_ProbChange_qualshock_all(w,1)
Av_Marg_Effect_qual_SE(w,1) = sqrt( Mean_ProbChange_qualshock_all_D2(w,1) * standard_errors(w,1)^2 * Mean_ProbChange_qualshock_all_D2(w,1) )

end

Price_sig = zeros(50,1);
Qual_sig  = zeros(50,1);
for w = 1:50
  if sign(Av_Marg_Effect_price(w,1) + 1.98*Av_Marg_Effect_price_SE(w,1)) == sign( Av_Marg_Effect_price(w,1))
        if sign(Av_Marg_Effect_price(w,1) - 1.98*Av_Marg_Effect_price_SE(w,1)) == sign( Av_Marg_Effect_price(w,1))
	  Price_sig(w,1) = 1;
	end
  end
  if sign(Av_Marg_Effect_qual(w,1) + 1.98*Av_Marg_Effect_qual_SE(w,1)) == sign( Av_Marg_Effect_qual(w,1))
        if sign(Av_Marg_Effect_qual(w,1) - 1.98*Av_Marg_Effect_qual_SE(w,1)) == sign( Av_Marg_Effect_qual(w,1))
	  Qual_sig(w,1) = 1;
	end
  end
end


[Av_Marg_Effect_price Av_Marg_Effect_price_SE Price_sig Av_Marg_Effect_qual  Av_Marg_Effect_qual_SE Qual_sig]
[Av_Marg_Effect_price Av_Marg_Effect_price_SE Av_Marg_Effect_qual  Av_Marg_Effect_qual_SE]

sum([Price_sig Qual_sig])

load ../china/data/trade_share.raw; 
trade_share=trade_share/sum(trade_share);

mean_Av_Marg_Effect_price=mean(Av_Marg_Effect_price);
Av_Marg_Effect_price_weighted=Av_Marg_Effect_price.*trade_share;
Av_Marg_Effect_price_weighted=sum(Av_Marg_Effect_price_weighted,1);
Av_Marg_Effect_price_weighted
mean_Av_Marg_Effect_price

mean_Av_Marg_Effect_price_SE=mean(Av_Marg_Effect_price_SE);
Av_Marg_Effect_price_SE_weighted=Av_Marg_Effect_price_SE.*trade_share;
Av_Marg_Effect_price_SE_weighted=sum(Av_Marg_Effect_price_SE_weighted,1);
Av_Marg_Effect_price_SE_weighted
mean_Av_Marg_Effect_price_SE

mean_Av_Marg_Effect_qual=mean(Av_Marg_Effect_qual);
Av_Marg_Effect_qual_weighted=Av_Marg_Effect_qual.*trade_share;
Av_Marg_Effect_qual_weighted=sum(Av_Marg_Effect_qual_weighted,1);
Av_Marg_Effect_qual_weighted
mean_Av_Marg_Effect_qual

mean_Av_Marg_Effect_qual_SE=mean(Av_Marg_Effect_qual_SE);
Av_Marg_Effect_qual_SE_weighted=Av_Marg_Effect_qual_SE.*trade_share;
Av_Marg_Effect_qual_SE_weighted=sum(Av_Marg_Effect_qual_SE_weighted,1);
Av_Marg_Effect_qual_SE_weighted
mean_Av_Marg_Effect_qual_SE

cd ../china/disclosure
save('coefficient_size','Av_Marg_Effect_price', 'Av_Marg_Effect_price_SE', 'Price_sig', 'Av_Marg_Effect_qual', 'Av_Marg_Effect_qual_SE', 'Qual_sig')

coefficient_size=[Av_Marg_Effect_price Av_Marg_Effect_price_SE Av_Marg_Effect_qual Av_Marg_Effect_qual_SE]
xlswrite('coefficient_size',coefficient_size)